NGB project: data analysis of bacterial and fungal data on Q. ilex trees
Marine Cambon
26 février, 2020
Data loading
The file "~/Documents/Postdoc_Biogeco/1_NGB/data/ORPHEE_leaf_microbiota/data_Qi_Qr_16S_filtered_formatted.RDS" and "~/Documents/Postdoc_Biogeco/1_NGB/data/ORPHEE_leaf_microbiota/data_Qi_Qr_ITS_filtered_formatted.RDS" are list objects with five elements:
samplesare the samples metadataasv_tableis a ASV x samples matrix, containing read counts of ASV per samplesassignis the taxonomic assignment of each ASVasv_sequencesis a fasta file with ASV sequences associated to their name in the other datasetswater_potentialis the predawn and midday water potential that have been measured concomitantly with leaf sampling for microbiota analysis
rm(list=ls())
# Load packages
library(vegan)
library(beanplot)
library(ape)
library(RColorBrewer)
library(knitr)
library(kableExtra)
library(ranger)
source("~/Documents/Postdoc_Biogeco/1_NGB/src/src_function_points_jitter.R")
source("~/Documents/Postdoc_Biogeco/2_DROUGHT/src/src_function_conv_hull_pcoa.R")
paf <- "~/Documents/Postdoc_Biogeco/1_NGB"
## Loading 16S data -------------
tmp <- readRDS("~/Documents/Postdoc_Biogeco/1_NGB/data/ORPHEE_leaf_microbiota/data_Qi_Qr_16S_filtered_formatted.RDS")
asv_table_16S <- tmp[["asv_table"]]
samples_16S <- tmp[["samples"]]
assign_16S <- tmp[["assign"]]
wp_16S <- tmp[["water_potential"]]
rm(tmp)
# Remove Qr samples from the dataset:
asv_table_16S <- asv_table_16S[,-grep("PU2", colnames(asv_table_16S))]
# Remove ASVs not present in samples:
asv_table_16S <- asv_table_16S[apply(asv_table_16S,1,sum)>0,]
# Remove assignation of ASVs not present in samples:
assign_16S <- assign_16S[match(rownames(asv_table_16S),
rownames(assign_16S)),]
# Remove Qr samples from the metadata
samples_16S <- samples_16S[match(colnames(asv_table_16S),
rownames(samples_16S)),]
# remove endophytes samples from the dataset:
endo_16S <- subset(samples_16S, sample_type=="endo")
asv_table_endo_16S <- asv_table_16S[,match(rownames(endo_16S),
dimnames(asv_table_16S)[[2]])]
samples_16S <- subset(samples_16S, sample_type=="micro")
asv_table_16S <- asv_table_16S[,match(rownames(samples_16S),
dimnames(asv_table_16S)[[2]])]
## Loading ITS data -------------
tmp <- readRDS("~/Documents/Postdoc_Biogeco/1_NGB/data/ORPHEE_leaf_microbiota/data_Qi_Qr_ITS_filtered_formatted.RDS")
asv_table_ITS <- tmp[["asv_table"]]
samples_ITS <- tmp[["samples"]]
assign_ITS <- tmp[["assign"]]
wp_ITS <- tmp[["water_potential"]]
rm(tmp)
# Remove Qr samples from the dataset:
asv_table_ITS <- asv_table_ITS[,-grep("PU2", colnames(asv_table_ITS))]
# Remove ASVs not present in samples:
asv_table_ITS <- asv_table_ITS[apply(asv_table_ITS,1,sum)>0,]
# Remove assignation of ASVs not present in samples:
assign_ITS <- assign_ITS[match(rownames(asv_table_ITS),
rownames(assign_ITS)),]
# Remove Qr samples from the metadata
samples_ITS <- samples_ITS[match(colnames(asv_table_ITS),
rownames(samples_ITS)),]
# remove endophytes samples from the dataset:
endo_ITS <- subset(samples_ITS, sample_type=="endo")
asv_table_endo_ITS <- asv_table_ITS[,match(rownames(endo_ITS),
dimnames(asv_table_ITS)[[2]])]
samples_ITS <- subset(samples_ITS, sample_type=="micro")
asv_table_ITS <- asv_table_ITS[,match(rownames(samples_ITS),
dimnames(asv_table_ITS)[[2]])]
# Colors
source("~/Documents/Postdoc_Biogeco/1_NGB/src/project_colors_v2.R")Sample metadata are loaded twice because samples that were successfully sequenced and that passed data curation are not exactly the same for 16S and ITS data. For now I don’t want to reduce the data sets to their intersection so it is more convenient to keep these duplicated data, but I could think later of a way to merge them.
Alpha diversity
Calculation of classical alpha diversity indices for each sample
samples_16S$shannon <- diversity(asv_table_16S, index = "shannon", MARGIN = 2)
samples_16S$chao1 <- estimateR(t(asv_table_16S))["S.chao1",]
samples_ITS$shannon <- diversity(asv_table_ITS, index = "shannon", MARGIN = 2)
samples_ITS$chao1 <- estimateR(t(asv_table_ITS))["S.chao1",]par(mfrow=c(1,2), bty="l",las=1, mar=c(4,4.5,2,0))
beanplot(chao1~irrigation,
data=samples_16S,
log="",
xlab="", col.axis="darkgray",
what=c(0,1,1,0),
cex.lab=1.5,
cutmin=0,
cex.axis=1,
ylim=c(0,300),
border="gray",
#border=colors[c("Bp", "Pp", "Qi", "Qr")],
#names=c("BRAN", "CDS1","CDS2","CDS3", "MTF1", "MTF2","MTF3"),
ylab = "Shannon index",
col=as.list("white"),
main="Bacteria")## new functionpt.jitter(samples_16S$irrigation, samples_16S$chao1,
bg = hue_irr[samples_16S$irrigation], alpha = 0.6)
beanplot(chao1~irrigation,
data=samples_ITS,
log="",
xlab="", col.axis="darkgray",
what=c(0,1,1,0),
cex.lab=1.5,
cutmin=0,
cex.axis=1,
ylim=c(0,300),
border="gray",
#border=colors[c("Bp", "Pp", "Qi", "Qr")],
#names=c("BRAN", "CDS1","CDS2","CDS3", "MTF1", "MTF2","MTF3"),
ylab = "Shannon index",
col=as.list("white"),
main="Fungi")## new functionpt.jitter(samples_ITS$irrigation, samples_ITS$chao1,
bg = hue_irr[samples_ITS$irrigation], alpha = 0.6)
par(mfrow=c(1,2), bty="l",las=1, mar=c(4,4.5,2,0))
beanplot(shannon~irrigation,
data=samples_16S,
log="",
xlab="", col.axis="darkgray",
what=c(0,1,1,0),
cex.lab=1.5,
cutmin=0,
cex.axis=1,
ylim=c(0,6),
border="gray",
#border=colors[c("Bp", "Pp", "Qi", "Qr")],
#names=c("BRAN", "CDS1","CDS2","CDS3", "MTF1", "MTF2","MTF3"),
ylab = "Shannon index",
col=as.list("white"),
main="Bacteria")## new functionpt.jitter(samples_16S$irrigation, samples_16S$shannon,
bg = hue_irr[samples_16S$irrigation], alpha = 0.6)
beanplot(shannon~irrigation,
data=samples_ITS,
log="",
xlab="", col.axis="darkgray",
what=c(0,1,1,0),
cex.lab=1.5,
cutmin=0,
cex.axis=1,
ylim=c(0,6),
border="gray",
#border=colors[c("Bp", "Pp", "Qi", "Qr")],
#names=c("BRAN", "CDS1","CDS2","CDS3", "MTF1", "MTF2","MTF3"),
ylab = "Shannon index",
col=as.list("white"),
main="Fungi")## new functionpt.jitter(samples_ITS$irrigation, samples_ITS$shannon,
bg = hue_irr[samples_ITS$irrigation], alpha = 0.6)
par(mfrow=c(1,2), bty="l",las=1, mar=c(4,4.5,2,0),
cex.lab=1.5, col.axis="darkgray")
plot(chao1~seq_depth,
col="gray",
pch=21,
bg=adjustcolor("gray", alpha.f = 0.6),
main="Bacteria",
data=samples_16S)
plot(chao1~seq_depth,
col="gray",
pch=21,
bg=adjustcolor("gray", alpha.f = 0.6),
main="Fungi",
data=samples_ITS)
plot(shannon~seq_depth,
col="gray",
pch=21,
bg=adjustcolor("gray", alpha.f = 0.6),
main="Bacteria",
data=samples_16S)
plot(shannon~seq_depth,
col="gray",
pch=21,
bg=adjustcolor("gray", alpha.f = 0.6),
main="Fungi",
data=samples_ITS)
The chao1 index clearly depends on sequencing depth. This should be taken into account somehow.
A sample has an extremely high number of reads in ITS run. I am not sure about what to do with it.
Community composition
Complete PERMANOVA
PERMANOVA analysis based on Bray-Curtis distance.
Bacteria
z <- samples_16S
tab <- asv_table_16S
(p <- adonis2(t(tab)~seq_depth+irrigation/block*month*leaf_age*psi_predawn_mean_sampling,
permutations = 999, method="bray",
data=z))
saveRDS(p, "PERMANOVA_Qi_tot_16S.RDS")
cat(
kable(p[,c(1,3,5)], digits = c(0,2,3),"latex", booktabs = T)%>%
row_spec(c(1:6), bold = T, color = "black"),
file ="tab_PERMANOVA_Qi_tot_16S.txt"
)readRDS("PERMANOVA_Qi_tot_16S.RDS")## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = t(tab) ~ seq_depth + irrigation/block * month * leaf_age * psi_predawn_mean_sampling, data = z, permutations = 999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## seq_depth 1 4.521 0.06200 17.1597 0.001 ***
## irrigation 1 1.964 0.02694 7.4551 0.001 ***
## month 1 1.322 0.01813 5.0172 0.001 ***
## leaf_age 1 4.118 0.05647 15.6286 0.001 ***
## psi_predawn_mean_sampling 1 0.637 0.00874 2.4177 0.005 **
## irrigation:block 6 2.543 0.03487 1.6085 0.001 ***
## irrigation:month 1 0.384 0.00526 1.4564 0.080 .
## irrigation:leaf_age 1 0.402 0.00551 1.5254 0.055 .
## irrigation:psi_predawn_mean_sampling 1 0.265 0.00363 1.0044 0.424
## month:psi_predawn_mean_sampling 1 0.208 0.00286 0.7909 0.727
## leaf_age:psi_predawn_mean_sampling 1 0.306 0.00419 1.1610 0.227
## irrigation:month:block 1 0.176 0.00242 0.6691 0.894
## irrigation:leaf_age:block 3 0.743 0.01019 0.9398 0.585
## Residual 210 55.333 0.75879
## Total 230 72.923 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Irrigation treatment explains 3% of the variance in bacterial community composition
- There is a significant effect of predawn water potential, but as we only have one value of wp per block, we need to disentangle the block and the irrigation effect. We should probably replace the block variable per a pcnm of the trees coordinates.
Fungi
z <- samples_ITS
tab <- asv_table_ITS
(p <- adonis2(t(tab)~seq_depth+irrigation/block*month*leaf_age*psi_predawn_mean_sampling,
permutations = 999, method="bray",
data=z))
saveRDS(p, "PERMANOVA_Qi_tot_ITS.RDS")
cat(
kable(p[,c(1,3,5)], digits = c(0,2,3),"latex", booktabs = T)%>%
row_spec(c(1:8, 11:13), bold = T, color = "black"),
file ="tab_PERMANOVA_Qi_tot_ITS.txt"
)readRDS("PERMANOVA_Qi_tot_ITS.RDS")## Permutation test for adonis under reduced model
## Terms added sequentially (first to last)
## Permutation: free
## Number of permutations: 999
##
## adonis2(formula = t(tab) ~ seq_depth + irrigation/block * month * leaf_age * psi_predawn_mean_sampling, data = z, permutations = 999, method = "bray")
## Df SumOfSqs R2 F Pr(>F)
## seq_depth 1 1.960 0.02510 7.6098 0.001 ***
## irrigation 1 1.619 0.02073 6.2857 0.001 ***
## month 1 3.086 0.03951 11.9790 0.001 ***
## leaf_age 1 7.268 0.09305 28.2124 0.001 ***
## psi_predawn_mean_sampling 1 0.670 0.00858 2.6010 0.002 **
## irrigation:block 6 5.121 0.06556 3.3129 0.001 ***
## irrigation:month 1 0.443 0.00567 1.7193 0.039 *
## irrigation:leaf_age 1 0.635 0.00813 2.4663 0.006 **
## irrigation:psi_predawn_mean_sampling 1 0.259 0.00332 1.0063 0.418
## month:psi_predawn_mean_sampling 1 0.278 0.00356 1.0780 0.320
## leaf_age:psi_predawn_mean_sampling 1 0.644 0.00825 2.5000 0.006 **
## irrigation:month:block 1 0.485 0.00621 1.8839 0.020 *
## irrigation:leaf_age:block 3 1.282 0.01641 1.6585 0.008 **
## Residual 211 54.356 0.69592
## Total 231 78.107 1.00000
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1- Irrigation treatment explains 3% of the variance in fungal community composition
- There is a significant effect of predawn water potential, but as we only have one value of wp per block, we need to disentangle the block and the irrigation effect. We should probably replace the block variable per a pcnm of the trees coordinates.
- It looks that there are more spatial effects in fungi than in bacteria
Visualization of the irrigation treatment effect
Clr transformation of the data
Due to the compositional nature of microbiota data and differences in sequencing depth between species, it is necessary to normalize data. I perform here a centered log-ratio (clr) transformation after modeling low count abundances with Monte Carlo replicates of the Dirichlet distribution (see Fernandes et al. (2014) Microbiome). The clr transformation shifts discrete data to continuous data.
Here 128 Monte Carlo samples of the Dirichlet distribution were made for each samples. The 128 matrices obtained were then averaged, following C. Pauvert’s thesis.
library(ALDEx2)
asv.clr <- aldex.clr(asv_table_16S, verbose=T,
useMC = F, denom = "all", mc.samples = 128)
asv.clr.mean<-sapply(getMonteCarloInstances(asv.clr),rowMeans) # Average over the different instances
saveRDS(asv.clr.mean,"aldex_clr_mean_16S.RDS") # Write to disk
asv.clr <- aldex.clr(asv_table_ITS, verbose=T,
useMC = F, denom = "all", mc.samples = 128)
asv.clr.mean<-sapply(getMonteCarloInstances(asv.clr),rowMeans) # Average over the different instances
saveRDS(asv.clr.mean,"aldex_clr_mean_ITS.RDS") # Write to diskBacteria
par(mfrow=c(1,2))
z <- samples_16S
tab <- asv_table_16S
dist <- vegdist(decostand(t(tab),method="hell"),
method = "bray", binary=F)
pp <- pcoa(dist)
coul <- adjustcolor(hue_irr[z$irrigation], alpha.f = 0.6)
# pdf(file.path(paf, "results/figures/fig_PCoA_hell_bray_Paracou_organ_16S.pdf"), 5, 5)
par(bty="l",las=1, mar=c(4,4,2,1))
plot(pp$vectors[,1], pp$vectors[,2], col=coul,
bg=coul,
pch=21, cex=1.5,
col.axis="grey50",
cex.lab=1,
cex.axis=1,
xlab=paste("PCoA axis 1 (",round(pp$values[1,"Relative_eig"],digits = 2)*100," %)", sep=""),
ylab = paste("PCoA axis 2 (",round(pp$values[2,"Relative_eig"],digits = 2)*100," %)", sep=""),
main="Hellinger, Bray-Curtis")
legend("bottomleft", pch=21,
pt.bg = hue_irr,
col=hue_irr,
legend=names(hue_irr))
conv.hull.pcoa(pp, groups = z$irrigation, colors=hue_irr)
# dev.off()
z <- samples_16S
tab <- readRDS("aldex_clr_mean_16S.RDS")
dist <- vegdist(t(tab), method = "euclidian", binary=F)
pp <- pcoa(dist)
coul <- adjustcolor(hue_irr[z$irrigation], alpha.f = 0.6)
# pdf(file.path(paf, "results/figures/fig_PCoA_hell_bray_Paracou_organ_16S.pdf"), 5, 5)
par(bty="l",las=1, mar=c(4,4,2,1))
plot(pp$vectors[,1], pp$vectors[,2], col=coul,
bg=coul,
pch=21, cex=1.5,
col.axis="grey50",
cex.lab=1,
cex.axis=1,
xlab=paste("PCoA axis 1 (",round(pp$values[1,"Relative_eig"],digits = 2)*100," %)", sep=""),
ylab = paste("PCoA axis 2 (",round(pp$values[2,"Relative_eig"],digits = 2)*100," %)", sep=""),
main="Clr, Euclidian")
legend("bottomleft", pch=21,
pt.bg = hue_irr,
col=hue_irr,
legend=names(hue_irr))
conv.hull.pcoa(pp, groups = z$irrigation, colors=hue_irr)
# dev.off()par(mfrow=c(1,2))
z <- samples_16S
tab <- asv_table_16S
dist <- vegdist(decostand(t(tab),method="hell"),
method = "bray", binary=F)
pp <- pcoa(dist)
coul <- adjustcolor(hue_leaf_age[z$leaf_age], alpha.f = 0.6)
# pdf(file.path(paf, "results/figures/fig_PCoA_hell_bray_Paracou_organ_16S.pdf"), 5, 5)
par(bty="l",las=1, mar=c(4,4,2,1))
plot(pp$vectors[,1], pp$vectors[,2], col=coul,
bg=coul,
pch=21, cex=1.5,
col.axis="grey50",
cex.lab=1,
cex.axis=1,
xlab=paste("PCoA axis 1 (",round(pp$values[1,"Relative_eig"],digits = 2)*100," %)", sep=""),
ylab = paste("PCoA axis 2 (",round(pp$values[2,"Relative_eig"],digits = 2)*100," %)", sep=""),
main="Hellinger, Bray-Curtis")
legend("bottomleft", pch=21,
pt.bg = hue_leaf_age,
col=hue_leaf_age,
legend=names(hue_leaf_age))
conv.hull.pcoa(pp, groups = z$leaf_age, colors=hue_leaf_age)
# dev.off()
z <- samples_16S
tab <- readRDS("aldex_clr_mean_16S.RDS")
dist <- vegdist(t(tab), method = "euclidian", binary=F)
pp <- pcoa(dist)
coul <- adjustcolor(hue_leaf_age[z$leaf_age], alpha.f = 0.6)
# pdf(file.path(paf, "results/figures/fig_PCoA_hell_bray_Paracou_organ_16S.pdf"), 5, 5)
par(bty="l",las=1, mar=c(4,4,2,1))
plot(pp$vectors[,1], pp$vectors[,2], col=coul,
bg=coul,
pch=21, cex=1.5,
col.axis="grey50",
cex.lab=1,
cex.axis=1,
xlab=paste("PCoA axis 1 (",round(pp$values[1,"Relative_eig"],digits = 2)*100," %)", sep=""),
ylab = paste("PCoA axis 2 (",round(pp$values[2,"Relative_eig"],digits = 2)*100," %)", sep=""),
main="Hellinger, Bray-Curtis")
legend("bottomleft", pch=21,
pt.bg = hue_leaf_age,
col=hue_leaf_age,
legend=names(hue_leaf_age))
conv.hull.pcoa(pp, groups = z$leaf_age, colors=hue_leaf_age)
# dev.off()par(mfrow=c(1,2))
z <- samples_16S
tab <- asv_table_16S
dist <- vegdist(decostand(t(tab),method="hell"),
method = "bray", binary=F)
pp <- pcoa(dist)
coul <- adjustcolor(hue_leaf_type[z$leaf_type], alpha.f = 0.6)
# pdf(file.path(paf, "results/figures/fig_PCoA_hell_bray_Paracou_organ_16S.pdf"), 5, 5)
par(bty="l",las=1, mar=c(4,4,2,1))
plot(pp$vectors[,1], pp$vectors[,2], col=coul,
bg=coul,
pch=21, cex=1.5,
col.axis="grey50",
cex.lab=1,
cex.axis=1,
xlab=paste("PCoA axis 1 (",round(pp$values[1,"Relative_eig"],digits = 2)*100," %)", sep=""),
ylab = paste("PCoA axis 2 (",round(pp$values[2,"Relative_eig"],digits = 2)*100," %)", sep=""),
main="Hellinger, Bray-Curtis")
legend("bottomleft", pch=21,
pt.bg = hue_leaf_type,
col=hue_leaf_type,
legend=names(hue_leaf_type))
#conv.hull.pcoa(pp, groups = z$leaf_type, colors=hue_leaf_type)
# dev.off()
z <- samples_16S
tab <- readRDS("aldex_clr_mean_16S.RDS")
dist <- vegdist(t(tab), method = "euclidian", binary=F)
pp <- pcoa(dist)
coul <- adjustcolor(hue_leaf_type[z$leaf_type], alpha.f = 0.6)
# pdf(file.path(paf, "results/figures/fig_PCoA_hell_bray_Paracou_organ_16S.pdf"), 5, 5)
par(bty="l",las=1, mar=c(4,4,2,1))
plot(pp$vectors[,1], pp$vectors[,2], col=coul,
bg=coul,
pch=21, cex=1.5,
col.axis="grey50",
cex.lab=1,
cex.axis=1,
xlab=paste("PCoA axis 1 (",round(pp$values[1,"Relative_eig"],digits = 2)*100," %)", sep=""),
ylab = paste("PCoA axis 2 (",round(pp$values[2,"Relative_eig"],digits = 2)*100," %)", sep=""),
main="Hellinger, Bray-Curtis")
legend("bottomleft", pch=21,
pt.bg = hue_leaf_type,
col=hue_leaf_type,
legend=names(hue_leaf_type))
#conv.hull.pcoa(pp, groups = z$leaf_type, colors=hue_leaf_type)
# dev.off()Fungi
par(mfrow=c(1,2))
z <- samples_ITS
tab <- asv_table_ITS
dist <- vegdist(decostand(t(tab),method="hell"),
method = "bray", binary=F)
pp <- pcoa(dist)
coul <- adjustcolor(hue_irr[z$irrigation], alpha.f = 0.6)
# pdf(file.path(paf, "results/figures/fig_PCoA_hell_bray_Paracou_organ_ITS.pdf"), 5, 5)
par(bty="l",las=1, mar=c(4,4,2,1))
plot(pp$vectors[,1], pp$vectors[,2], col=coul,
bg=coul,
pch=21, cex=1.5,
col.axis="grey50",
cex.lab=1,
cex.axis=1,
xlab=paste("PCoA axis 1 (",round(pp$values[1,"Relative_eig"],digits = 2)*100," %)", sep=""),
ylab = paste("PCoA axis 2 (",round(pp$values[2,"Relative_eig"],digits = 2)*100," %)", sep=""),
main="Hellinger, Bray-Curtis")
legend("bottomleft", pch=21,
pt.bg = hue_irr,
col=hue_irr,
legend=names(hue_irr))
conv.hull.pcoa(pp, groups = z$irrigation, colors=hue_irr)
# dev.off()
z <- samples_ITS
tab <- readRDS("aldex_clr_mean_ITS.RDS")
dist <- vegdist(t(tab), method = "euclidian", binary=F)
pp <- pcoa(dist)
coul <- adjustcolor(hue_irr[z$irrigation], alpha.f = 0.6)
# pdf(file.path(paf, "results/figures/fig_PCoA_hell_bray_Paracou_organ_ITS.pdf"), 5, 5)
par(bty="l",las=1, mar=c(4,4,2,1))
plot(pp$vectors[,1], pp$vectors[,2], col=coul,
bg=coul,
pch=21, cex=1.5,
col.axis="grey50",
cex.lab=1,
cex.axis=1,
xlab=paste("PCoA axis 1 (",round(pp$values[1,"Relative_eig"],digits = 2)*100," %)", sep=""),
ylab = paste("PCoA axis 2 (",round(pp$values[2,"Relative_eig"],digits = 2)*100," %)", sep=""),
main="Clr, Euclidian")
legend("bottomleft", pch=21,
pt.bg = hue_irr,
col=hue_irr,
legend=names(hue_irr))
conv.hull.pcoa(pp, groups = z$irrigation, colors=hue_irr)
# dev.off()par(mfrow=c(1,2))
z <- samples_ITS
tab <- asv_table_ITS
dist <- vegdist(decostand(t(tab),method="hell"),
method = "bray", binary=F)
pp <- pcoa(dist)
coul <- adjustcolor(hue_leaf_age[z$leaf_age], alpha.f = 0.6)
# pdf(file.path(paf, "results/figures/fig_PCoA_hell_bray_Paracou_organ_16S.pdf"), 5, 5)
par(bty="l",las=1, mar=c(4,4,2,1))
plot(pp$vectors[,1], pp$vectors[,2], col=coul,
bg=coul,
pch=21, cex=1.5,
col.axis="grey50",
cex.lab=1,
cex.axis=1,
xlab=paste("PCoA axis 1 (",round(pp$values[1,"Relative_eig"],digits = 2)*100," %)", sep=""),
ylab = paste("PCoA axis 2 (",round(pp$values[2,"Relative_eig"],digits = 2)*100," %)", sep=""),
main="Hellinger, Bray-Curtis")
legend("bottomleft", pch=21,
pt.bg = hue_leaf_age,
col=hue_leaf_age,
legend=names(hue_leaf_age))
conv.hull.pcoa(pp, groups = z$leaf_age, colors=hue_leaf_age)
# dev.off()
z <- samples_ITS
tab <- readRDS("aldex_clr_mean_ITS.RDS")
dist <- vegdist(t(tab), method = "euclidian", binary=F)
pp <- pcoa(dist)
coul <- adjustcolor(hue_leaf_age[z$leaf_age], alpha.f = 0.6)
# pdf(file.path(paf, "results/figures/fig_PCoA_hell_bray_Paracou_organ_16S.pdf"), 5, 5)
par(bty="l",las=1, mar=c(4,4,2,1))
plot(pp$vectors[,1], pp$vectors[,2], col=coul,
bg=coul,
pch=21, cex=1.5,
col.axis="grey50",
cex.lab=1,
cex.axis=1,
xlab=paste("PCoA axis 1 (",round(pp$values[1,"Relative_eig"],digits = 2)*100," %)", sep=""),
ylab = paste("PCoA axis 2 (",round(pp$values[2,"Relative_eig"],digits = 2)*100," %)", sep=""),
main="Hellinger, Bray-Curtis")
legend("bottomleft", pch=21,
pt.bg = hue_leaf_age,
col=hue_leaf_age,
legend=names(hue_leaf_age))
conv.hull.pcoa(pp, groups = z$leaf_age, colors=hue_leaf_age)
# dev.off()par(mfrow=c(1,2))
z <- samples_ITS
tab <- asv_table_ITS
dist <- vegdist(decostand(t(tab),method="hell"),
method = "bray", binary=F)
pp <- pcoa(dist)
coul <- adjustcolor(hue_leaf_type[z$leaf_type], alpha.f = 0.6)
# pdf(file.path(paf, "results/figures/fig_PCoA_hell_bray_Paracou_organ_16S.pdf"), 5, 5)
par(bty="l",las=1, mar=c(4,4,2,1))
plot(pp$vectors[,1], pp$vectors[,2], col=coul,
bg=coul,
pch=21, cex=1.5,
col.axis="grey50",
cex.lab=1,
cex.axis=1,
xlab=paste("PCoA axis 1 (",round(pp$values[1,"Relative_eig"],digits = 2)*100," %)", sep=""),
ylab = paste("PCoA axis 2 (",round(pp$values[2,"Relative_eig"],digits = 2)*100," %)", sep=""),
main="Hellinger, Bray-Curtis")
legend("bottomleft", pch=21,
pt.bg = hue_leaf_type,
col=hue_leaf_type,
legend=names(hue_leaf_type))
conv.hull.pcoa(pp, groups = z$leaf_type, colors=hue_leaf_type)
# dev.off()
z <- samples_ITS
tab <- readRDS("aldex_clr_mean_ITS.RDS")
dist <- vegdist(t(tab), method = "euclidian", binary=F)
pp <- pcoa(dist)
coul <- adjustcolor(hue_leaf_type[z$leaf_type], alpha.f = 0.6)
# pdf(file.path(paf, "results/figures/fig_PCoA_hell_bray_Paracou_organ_16S.pdf"), 5, 5)
par(bty="l",las=1, mar=c(4,4,2,1))
plot(pp$vectors[,1], pp$vectors[,2], col=coul,
bg=coul,
pch=21, cex=1.5,
col.axis="grey50",
cex.lab=1,
cex.axis=1,
xlab=paste("PCoA axis 1 (",round(pp$values[1,"Relative_eig"],digits = 2)*100," %)", sep=""),
ylab = paste("PCoA axis 2 (",round(pp$values[2,"Relative_eig"],digits = 2)*100," %)", sep=""),
main="Hellinger, Bray-Curtis")
legend("bottomleft", pch=21,
pt.bg = hue_leaf_type,
col=hue_leaf_type,
legend=names(hue_leaf_type))
conv.hull.pcoa(pp, groups = z$leaf_type, colors=hue_leaf_type)
# dev.off()Session info
sessionInfo()## R version 3.6.1 (2019-07-05)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Debian GNU/Linux 9 (stretch)
##
## Matrix products: default
## BLAS: /usr/lib/libblas/libblas.so.3.7.0
## LAPACK: /usr/lib/lapack/liblapack.so.3.7.0
##
## locale:
## [1] LC_CTYPE=fr_FR.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=fr_FR.UTF-8 LC_COLLATE=fr_FR.UTF-8
## [5] LC_MONETARY=fr_FR.UTF-8 LC_MESSAGES=fr_FR.UTF-8
## [7] LC_PAPER=fr_FR.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ranger_0.12.1 kableExtra_1.1.0 knitr_1.28 RColorBrewer_1.1-2
## [5] ape_5.3 beanplot_1.2 vegan_2.5-6 lattice_0.20-38
## [9] permute_0.9-5
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.3 compiler_3.6.1 pillar_1.4.3 tools_3.6.1
## [5] digest_0.6.24 viridisLite_0.3.0 lifecycle_0.1.0 evaluate_0.14
## [9] tibble_2.1.3 nlme_3.1-140 mgcv_1.8-28 pkgconfig_2.0.3
## [13] rlang_0.4.4 Matrix_1.2-17 rstudioapi_0.11 yaml_2.2.1
## [17] parallel_3.6.1 xfun_0.12 xml2_1.2.2 stringr_1.4.0
## [21] httr_1.4.1 cluster_2.1.0 vctrs_0.2.2 hms_0.5.3
## [25] webshot_0.5.2 grid_3.6.1 glue_1.3.1 R6_2.4.1
## [29] rmarkdown_2.1 readr_1.3.1 magrittr_1.5 scales_1.1.0
## [33] htmltools_0.4.0 MASS_7.3-51.4 splines_3.6.1 rvest_0.3.5
## [37] colorspace_1.4-1 stringi_1.4.5 munsell_0.5.0 crayon_1.3.4